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ABSTRACT 


A three dimensional staggered grid model, based on the Boussinesq 
equations, is developed, and applied to mesoscale waves. The model is 
tested with a three dimensional stable gravity wave and forecasts for 
40 time steps are compared with the analytical solutions and also with 
a full grid model. Comparisons show very good agreement with both 
forecasts for the stated time. Failure of the relaxation scheme to 


meet the convergence criterion after 47 time steps is discussed. 
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IT. INTRODUCTION 


In the past few years many attempts have been made to understand 
the complex dynamic interactions between the thunderstorm and its 
environment. Newton and Newton (1959) proposed a dynamical thunder- 
storm model which proposed that strong vertical motions within the 
thunderstorm cell create an effective barrier to the flow of the 
environmental air. 

Ogura (1963) first discussed the probability of using the primitive 
equations for forecasting small scale dynamical features. Interaction 
between the convective elements and the environmental wind field can 
be investigated numerically by the use of these equations. Deardorff 
(1970) made a numerical study of turbulent flow using a three-dimensional 
primitive equation model. A similar forecasting procedure is used in 
this model. The importance of gravity waves in small scale convective 
activity is discussed by Matsumoto and Ninomiya (1969). Orszag (1969) 
presented a model for the simulation of turbulence which employs a 
fast Fourier Transform instead of the Liebmann Relaxation Technique. 
Further refinement of this model may be accomplished by the use of 
such a method. 

Norton (1970) discussed the development of a three dimensional model 
for application to mesoscale motions based on the Boussinesq equations. 
The model presented here is a further development of the same model with 
a staggered grid system applied to minimize computer storage requirements. 
The total grid size is 14 by 14 in the horizontal with 10 levels in the 


vertical, while the computational grid size is 11 by 10 by 6. The grid 





interval in all directions is |] kilometer. The Coriolis parameter, 

treated as a constant, is determined at a mean latitude of 35° North. 
Numerical computations are made using the IBM 0S/360 Computer of 

the W. R. Church Computer Center. To conform to computational stability 


a time step of 30 seconds was chosen. 





IT. MODEL DEVELOPMENT EQUATIONS 


A. FORECAST EQUATIONS 


The pressure term in the complete equation of motion can be written 


as: 
] RT Po 1/« 
p 
or 
lt ypzerv g 
Gres 3 
where 
K 
b= 6, = 
Pi Py 
and Z 
p 
a= 1 {79 
p 
Using these results in the component form of the equation of motion 
yields: 
duo. _ 9 38 Z 
at aN eR (al) 
dv . _ 9 38 2 
dt a a W We 3 (2) 
dv 4 8 | 2 
a 9 ay g tw4w . (3) 


The perturbation expressions can be formed by: 


8 = 6 (z) +8' (x,y,z,t) 


0 
6 = 8 + 8(xsysZ,t) 
where 
Ig'| << B and lo | Say 


and the latter inequality corresponds to the Boussinesq approximation. 
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After substitution into equations (1), (2), and (3), the equations 


of motion can be written: 


Gl) ee 06 2 
Ree (Y  (YS UD (4) 
dv _ _, 38° | 2 
dt °6 by ES AE: (5) 
dB dB 
dws w0: OI as eels ye 
and ae) 5 =a 3 5s g tvV¥3 W . 
with the relatively small products: 
ap’, at ap" 
0 aX 5 8 3y » and 6 = 
being neglected. The hydrostatic equation relates By and 8 by: 
2 ete 
JZ By 
and the vertical component can now be written: 
dw. _, 38' , 98 42 
dt ry i 3 sae! (6) 
By defining ¢ = 6,8' and substituting into (4), (5), and (6) and 
expanding the total derivatives, the system of equations becomes: 
cnn ih eee 2 
= LAUDS ees TE Sie Te (7) 
ae Se ae 2 
L(v) v7 fu + v¥3 Vv, (8) 
ON ae _ 99 go 2 
at LED) es G wae (9) 
96 | 
Uae V = 0 (ile) 
where L(S) = - V3° (V4S), S = U,v.w.6 . 


1] 





Incompressible flow is assumed, to eliminate sound waves. This system 
has become known as the anelastic equations and are quite accurate for 


shallow convection. (Ogura and Phillips, 1962) 


B. DETERMINATION OF THE PRESSURE TERMS 
Equations (7), (8), and (9) may be combined in vector form and 
written: 
aV, a a oe =: qe - a 
ar ae -(V. . v)V3 - Va¢ -f(k x V3) sa kK + V4 V2 (12) 
O 
By taking the three dimensional divergence and assuming the time rate of 


change to be zero, (12) becomes: 


(V4 ° V3) - 9 


-V3 [(V3-¥3)V3] - v§9 
(13) 


For the horizontal extent used in this model, the Coriolis parameter may 


be treated as a constant and then (13) becomes: 


2, = -y.¢[(Ve-v.)V (ie x) we 
eg 1GYGIV3) > 7, va°(k x V3) + a, 8 
— = e i on 
¥3°[(V3+73)V3] + fF, c + @ 32 +¥3 +VV3V (14) 


Note that when the terms on the RHS are solved by ordinary differencing 
techniques , the value of ¢ on the LHS can be solved by using the 


Extrapolated Liebmann Relaxation Scheme. 
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III. FINITE DIFFERENCING 


Numerical integration of equations (7) through (10) is carried out 
using finite difference methods in x,y,z, and t. Extrapolated Liebmann 
Relaxation is used to solve for the pressure term in equation (14). For 
the complete development of the relaxation scheme, see Appendix A. The 
technique applied to generate the forcing function used in (14) is shown 
in complete form in Appendix B. 

To maintain consistency in finite differencing throughout the fore- 
cast procedure, the differencing scheme used in equation (14) must be 
equivalent to the type used on equations (7), (8), and (9). The flux 
form differencing method devised by Arakawa (1966) was utilized in the 


form: 


BCS) 


is5sk © loss eS 


ef oso? in le 
SU aeguk detg sk 55544) f 4x 
+04 jerk PMs pe) (Ss gee § Saas 

a. oe elk * S534) | / A Ay 
los je + Wisk) (Sig gay + Si 5k) 

SWE gk + gk?) Sag ke1 * S530 | / Anz 


This method is used for the u,v,w and o fields, but for the linear space 
derivatives for 4, an averaging technique must be employed due to the 
Staggered grid system. The centered form of the space derivative takes 


on the form: 





9S - 
Fr), yo ; Is, »J sk y > k=! >: alok “ 4 el fell! / 4. 
1',J'»k 
tee eoaetgnke1) © ~i-1,5-1,k ° °t-1i,d~1,k-1//4-|/4% 


with similar forms for the y and z derivatives. The primed subscripts 
refer to the grid points one-half grid distance above or below the 

diagonals connecting surrounding velocity points in a horizontal plane. 
Centered time differencing ("leapfrog") is used for all forecasts with 


the exception of the first time step where a single forward time step 


is used. 


14 


— 





IV. BOUNDARY CONDITIONS 


The boundary conditions are selected to enforce mass and energy 
conservation within the confines of the computational grid. Figure | 
illustrates the horizontal grid domain and boundary placement. Cyclic 
continuity of parameters is imposed by forcing periodicity in the 


east-west direction according to the relations 


°1,.0,K ~ OIM-2,d,K 
309.K ~ OIM-1,0,K 
SIMJ.K ~ 53.0. 


where S = u,v,W, 6. Since values of S are required at three consecutive 
points in solving for the pressure and divergence terms, values outside 
the domain can be obtained through the cyclic conditions. The east- 


west boundary conditions on ¢ are also cyclic and are specified by: 


o10.K ~ ?IM-2,0,K 


OTN el [ee 


The south boundary is located halfway between velocity grid points 


J = 2 and J = 3, while the north boundary is between the points 
J = JM - 2 and J = JM- 1. To prevent the flow of mass and energy 
through these boundaries, the v-component is set by the following 


conditions: 





j=5M sa 
e e o e e ar JM 
4 a 6 4 
e e e | e e e JM-1 
North 
a 
9 9 r © r e JM-2 
A rN rN rN 
. e e r ® @ 
® @ © ® @ 6 
Aa a 4 a 
e e e © e 3 
South 
@ r) e e @ e 2 
Y a yal 4 & 
r ® e * © © 
i- 1 2 3 IM-2 {M-1 1= IM 
jr 1 a i= 1 
e-U V,W,8 


Figure 1. Illustration of the grid system in the horizontal 
direction. 
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V -V 


I,JM,K I,JM-3,K 


V -V 


I,JM-1,K I ,JM-2 ,K 


The values of u,w, and 6 are extrapolated outward by: 


SroakK 7 >1,2,K 7 °4,3.K 


S S S 


I,uM,K I,JM-1,K T,JM-2,K 


where S = u,w,6. This implies that there is no advection through this 
boundary and furthermore no heat exchange as discussed by Ukaji and 
Matsuno (1970). 

The geostrophic relationship: 


(5 


* 
= -f Us5 24 i 
ae Om lh aoa 


was used to obtain the values of » outside the north-south boundaries. 


* 
i'sj'ok' represents the volume average of the u velocity field at 


the boundary. The results reported here were for the special case of 


Where u 


f= 0. For this case, the boundary condition then becomes: 


O 
2) = 0 
5 i 


Figure 2 illustrates the vertical grid domain and boundary placements. 
The top and bottom boundaries are set in manner similar to the north- 
south boundaries where the points midway between K = KM - 1 and K = KM - 2 
and K = 2 and K = 3 are the top and bottom boundaries respectively. The 
no-flux condition is imposed on the w field by: 


W = -W 


st eg I,J,4 


W = -W 


les 52 AUS: 


17 








k=KM k=KM 
® ry 2 e © e 
PA rs S “ 
i a ie : e e KM-1 
Top 
4 
e e © @ € e KM-2 
mn n A PN 
© a 8 @ & e KM73 
4 A a r 
; . e 8 . * 
4 rN a & 
a e e e © * 
a A rs a 
e * e ; r @ e 3 
Bottom 


e @ a e ° * 2 
& A & A 
e € e e ° 
i 2 3 IM-2 1M -1 i— IM 
k= a-¢ k= 1 
°-U,V,W,2 


Figure 2. Illustration of the grid system in the vertical 
a direction. 
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W -W 


I,J ,KM I,J ,KM-3 


W -W 


1,J0,KM-1 I,J,KM-2  , 


while the 6 values are set by the relations: 


ae OR 
seen eee 1 .J,3 
OT d,KM 8a, KM-3 
ORK s 1  Ind, KM-2 


The values of u and v are simply extrapolated outward by: 


ee)  smeiedk? = 10.3 


S S 


1,3,KM Lali) © “yah 


where S = u,v. 
The value of the pressure term is obtained by the use of: 


Fae 
Sd dak 


The importance of proper boundary conditions cannot be over- 
emphasized. This model was found to be extremely sensitive to all the 
boundary conditions. The complete development of the boundary conditions 


1s found in Appendix C. 
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V. INITIAL CONDITIONS 


Initialization of the u,v,w and 6 fields 1s accomplished by the use 
of linearized gravity wave solutions, which are exact for small amplitude 


gravity waves and are derived from the following equations: 


ou. _ 8b 
at OX 
ove. _ 9¢ 
ee dy 
ow = = ao. + go 
ot 9Z ) 
0 
We id V = 0 
and _ 
90 | 30 
ot x OZ 


The exact solutions for a stable wave propagating in the x direction are 


given by: 

U =- as SMe =Ge)ecOS jy GOS Mz , 
2 2 
A ty 

vVo=- um cos(a MESChINSiINe py cos mz , 
2 a 
A + y 

w = Wcos(x x -ct) cos py sin mz , 

W 90 .. ; 
Os ane ay Sin(r x -ct) cos py sin mz 


where c represents the phase speed of the stable gravity wave in Mm sec 
and is calculated from: 


Soe pa ee (15) 
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and N is the Brunt Vaisala frequency 


The initial guess for the ¢ field is chosen to be an arbitrary constant 
(2.0) due to the extremely efficient relaxation technique. The set of 


exact solutions with t = 0 was used to initialize the fields. 
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Wie  REOUETS 


By introducing a stable gravity wave perturbation into the basic 
model, the wave form can be tested for proper phase speed, amplitude 
growth or decay and computational stability. The initial conditions 
should result in a wave that propagates in the x direction alone with 
no growth in amplitude. Figures 3a, 3b, 3c, and 3d are representative 
plotted values of the u,v,w, and 6 wave forms after time steps 10, 

20, 30, and 40, with the initial value present for comparison. These 
plots are for the wave forms in the x direction for J = 5 and K = 7 
and show the phase speed relationship between the various time steps. 


The exact solution to equation (15) yields a phase speed of 5.5 m sec. 


| ome le Ae Wer’ eae! 


] 


for the specified values W = .01 m sec. 
Phase speed calculated from the plotted values is 5.4 msec for each 

of the four wave forms. There is negligible amplitude growth or decay 

in any of the plots. 

Figures 4a, 4b, 4c, and 4d are the plotted values of the wave forms 
in the y direction for I = 5 and K = 7. The plots indicate there is 
no propagation in the y direction. Due to the nature of the solutions 
to the linearized equations, the wave forms should oscillate about the 
zero value as indicated in the stated figures. 

Plots of the wave forms in the z direction for I = 5 and J = 9 are 
Shown in Figures 5a, 5b, 5c, and 5d and are similar in amplitude 
oscillation to the y direction with no phase speed observed. 

The main objective of this staggered grid model was to reduce the 
storage size and computer time required for an operational small scale 


convective model. The previous full-grid model required 502K bytes of 
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Storage and approximately 20 minutes of computer time to produce a 40 
time step forecast. The forecasts produced by the staggered grid 
model accomplishes a 40 time step prognosis in less than 10 minutes of 
computer time and requires only 156K bytes of storage. 

The results obtained after 40 time steps compared to the analytical 
solutions illustrate the accuracy of the finite difference solutions. 
However, a few time steps later, the solutions fail to converge. A 


further discussion of this difficulty will be given in a later section. 
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VII. CONCLUSIONS AND SUMMARY 


A ten-level staggered grid primitive equation model was applied to a 
mesoscale grid and tested with an initial stable gravity wave perturba- 
tion to check the computational stability and finite difference form of 
the equations of motion. Comparing the results presented here to a non- 
staggered grid model, the computer time and storage requirements have 
been reduced to about 1/3 of the full grid model. 

The forecasts produced by this model compare quite closely with the 
full grid model but ceases to function after 47 time steps. Up until 
this time step all fields appear to be predicted properly. Tests were 
run in double precision to determine if truncation error was a problem 
and the same type of results occured with the model failing to converge 
after 47 time steps. In all cases the points that failed to converge 
were near the boundaries which leads to the conclusion that further 
refinements of the boundary conditions may be necessary. Experiments 
have shown that a less stringent relaxation tolerance produces less 
exact @ gradients which in turn induces greater errors in the prognosis 
of the other parameters. Apparently there is still some difficulty in 
the program and presently studies are being co. ducted to further refine 
this model. 

Continued studies using this model should have a provision for 
including vertical compressibility. Additionally, since a primary 
source of the energy in a thunderstorm-type cloud is latent heat, 
application to cumulus scale convective activity would be inadequate 


without the inclusion of moisture. 
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APPENDIX A 
EXTRAPOLATED LIEBMANN RELAXATION SCHEME 
FOR THREE DIMENSIONAL APPLICATION 


The equation to be solved by relaxation is: 


y@ -F = 0 (A-1) 


where the tendency of the divergence is zero and: 


_. id. ev gq. 20 ; 
Fos -V, (, 7473+ oe (A-2) 


with the ¢ points staggered in the vertical as well as in the horizontal. 
(They are located one-half grid distance above or below the diagonals 
connecting velocity points as in Fig. Al.) In other words, the velocity 
and 6 points are located at the corners of a cube and the ¢ point in the 
center. For convenience, a new temporary indexing scheme will be 


utilized by writing: 


i’ = i - u 
2 

1 

t = _ cee 

J J 9 
_ ae 
k' = k 5 


where i,j,k represent u,v,w and 6 points and i', j', k' represent 
points. Writing (A-1) with this new notation in finite difference form 
at the point i', j', k' and expanding yields: 
a + 2 + a. 
eee se tlejink sl) 741 ,5°-1,k'41 ~ PF", 547, kt 4 
(4ax) 
-4¢ 


Pee ey Tyg te ki4] © OFs'-1,9",kt41 
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Three dimensional schematic of the location of 


Various fields. 


Figure Al. 
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+ 26.1 + 4 


aU 9) bey) (ley t1,g't) kt * Ta'ag' kt * 9 '41,95'-1,k! 


ee ee Peg yk 41,9" 41,k' 


eek w ob 51g '-1,k! PFET tl kid * SO a Oe 


= 2 


Tee ete peg 8g! kK’) OP" sj '=1,k'=1 


eed ais ipkiot o°s'21,5",k'-1 * treagtetkit |? 


1 


Taye [tieta atk — bray a kia t Spar gteT klar 7 Og per ke 


Foe pt gtar F ebg pea kit FSET gered 7 MOET g kta 


+ el.) Swe a8 ob aT t+] k = O14] 4" k! + 2b 1aT j'-1,k! 


at UES Bee - Boar gt at wee SURG) Ai ce 2b 511 j' 41k! 


ae eM’ s sid lek) °t'41,g°41.k'-1 > “°4'41,5°,k'=1 


+ 265, 


TOPHL THT kt-1 PHT RT ORT K ETT EFT kta 


Meegiek eT) <Pi=1,g°,k°-1 Paget ke] . 


] 


(4az)° PP eT 41, k'+1 © ee aT gtk 4] : es eee oes ALES aes 


+4. 7 ob a +: ob 


ji sk'+ fteT kit) * P41, 5 '41 kt SS) 5.9 
sei tees aiey ike © ist, 5',k’ ~ °P4'41,9°-1,k' 
ee apie pte lak? o°4'-1,5'41,k' 


eee ese iedak) © fitei,j'4iyk'=1 © 2¢4'21,5",k'-1 


oy 





+ 26, 


Tete es y 41k 84", 5", ki see 


gene a) | CPa'=1,5",k'-1 * #91, 511 | 


- F 0 (A-3) 


eras 
Rearranging terms and clearing fractions leads to: 
P55 


+(4ax)* (4az)¢ c - iia 


(4ay)* (4az)é 


+(4ax)* (4ay)* c B05 5k 


-(4ax)© (day) (4az)* F 0 (A-4) 


ies sk 
where B, C and D represent the finite difference forms in equation (A-3) 
with the term Boa gtk! subtracted off. With further rearranging of 
terms €-4) becomes: 


Z Ze 2 


bar vi pi =A (aay)? (An Zee tides (44z)~ C 


TosJ > 


+ (4ax)* (day)? D - (4ax)? (day)® (4az)? Fanaut 


where: 
; 
A = 


2 (aay)? (4az)* + (4ax)® (40z)* + (4ax)? (aay)? | 
; old 
Adding and subtracting ba j'.k' yields: 


ld Z 2 
OFT gtk a Oph gh eaten tee) (4az)* B + (4ax)* (4az) Cs 


(A- 5) 


2 2 2 


+(4ax)” (4ay)" D - (4ax)? (4ay) 


2 
(4Az) ratake | 
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Where the over-relaxation coefficient a is introduced. (A-5) then 


becomes: 


new 7 old 2 2 
Oi 5k! oT (1 = a) oe ee 5 el 6 Ay (doy) (4az) B 


+ (4ax)® (4az)* C + (4ax)*® (4ay)® D 


- (Aax)® (4ay)*® (4az)* Finan | (A-6) 


This was the form used for relaxation in the model. Many tests were run 
to determine the most efficient relaxation coefficient for this type of 
finite difference scheme. The value of a which resulted in the most 

efficient relaxation was 1.05, Note that if Ax = Ay = Az, equation (A-6) 


reduces to the more familiar form: 


new _ old oO 
Op gt ge = (1 =o) ban gt Qe F 9g | B+c+D 
2 
7 (442) Fie gt og ° 


The stencil of the method of the ye operator in the three dimension- 
al staggered grid, assuming that the three grid jengths are equal, is 


shown in Fig. A2. 


4] 








¥ Operator in three dimensional form assuming 


Figure A2. 


AX = Ay = AZ. 
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APPENDIX B 
DEVELOPMENT OF THE FORCING FUNCTION 


The formation of the forcing function used in the relaxation scheme 
is a long process. Consistency between the finite difference methods 
used in the prognostic equations and in the relaxation technique must 
be maintained. The computation of the forcing function given in 
equation (A-2) in finite difference form follows. In order to determine 
this forcing function in a consistent manner, the forcing function must 
first be calculated at an i,j,k point and then volume averaged to 
have the function at an i',j',k' point. Each of the three directional 
Components are computed separately and then summed to put the complete 
forcing function at the desired point. For example, in order to form 
a Glee (the forcing function at the point i' = 5, j' = 5, k' = 5) 


the following steps are required: 


x J] 
“50 _0ntat “1 | 816, 6.6 + Glu), 5g + Glue go + Glude cg 


G(u) + G( 


RG AG We 56 * a 


Ws 6.5 


G(v) + G( 


S166 De Bg 


Demers 


Z Be: 
oes! 5° = 3 GlW6 6,6 * 26 ,5,6 * Swe gg * Sle c.g 


[e(v), 5,6 * Ve 5g * Ble 55 # a(V)5.5,5| 


alt 1 

ZG, 6 ot Gwe cet Gwe 6 et Glwie os Hi 
ES eX y z 

Fgh eit gt iceceteee Sciec iru? 8) ci Gils 
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where 


G(u) = - L(u) + fv t vw5u , 


G? 
~~ 
<= 
—- 
iT) 


- L(w) + 22 + vv3 w 
8 


G> 

~~ 
= 

—— 
i 


2 
3 
2 
- L(v) - fu + Ww Ve an 
2 
an 


and the G functions are calculated on the i,j,k grid. 
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APPENDIX C 
COMPLETE DEVELOPMENT OF THE BOUNDARY CONDITIONS 


For the boundary described by the plane of j' = 2° the y component 


of the equation of motion may be written 


oV 
ON = 0 
Obit Ook! 
which implies 
av - _ wv 
ati 34k a eek (C-1) 
Further expansion of each of the terms leads to: 
JV : ad 
= = L(y). - - fu. (C-2 
ate 3k es Ooi Wa k Omak 
and 
- se 2s am 2 as cet ae (C-3) 
ie ak ease ef 


By introducing the no-flux conditions on v at the lateral walls from 


(C-1), (C-2) and (C-3) yield: 


’ 


IS 


_ ou 


(C-4) 
°Y4 3,k 


- fu. = oe bie 
Omi. Kk aYs 2k Omdiiaek 


Expansion of (C-4) into finite differences and averaging the u 
field to be consistent with the > field yields: 


it 
in = 0 (u 


OF ask! > $4 (3k 1525k 7 Ui41,2,k 7 MG,2 KH 


pau ay 
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with the no-flux conditions being: 


a ao Kk 


V. ane 
Teleek =, Va 4,k 


A similar approach is taken at thenorthern boundary with the following 


conditions: 


viedaete CO ed ek 


ViUM.K 0 4,03,k 


_ _ Ay'o 
bear! = bGja3,k! 7 2 MY ig2yk * YG41,02,k 


Ue spend 


sf Napanee, 
The case considered here was for Ie = 0. 


For the boundary described by k' = 2', the z component of the 


equation of motion may be written: 


3W 
ow = 8 
ater on yt 
which implies: 
ow ee us (C-5) 
ati 3.k ati ok 


Further expansion of each of the terms in (C-5) yields: 


ow. See . 96 2 


0. (C-6) 
oti 3k Ie ace cone 12-2 
_ ow : 9d _ : 
at es te oe a a C-7) 
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By introducing the no-flux conditions on w at the top and bottom 


boundaries from (C-5), (C-6) and (C-7) they yield: 


(C-8) 


By setting © at the boundaries as previously discussed, (C-8) is 


satisfied if: 


At | en oe 
where the 6 boundary conditions are: 

ene aa 3 

gel) | aga 


A similar approach is taken at the top boundary with the following 


conditions: 

elas sis aK? 

PF pKM 785 ,5,K3 
and 

OF 5" KT = Ot 5! KB! 
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